Thermopower of Single- Molecule Devices 
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We investigate the thermopower of single molecules weakly coupled to metallic leads. We model 
the molecule in terms of the relevant electronic orbitals coupled to phonons corresponding to both 
internal vibrations and to oscillations of the molecule as a whole. The thermopower is computed 
by means of rate equations including both sequential-tunneling and cotunneling processes. Under 
certain conditions, the thermopower allows one to access the electronic and phononic excitation 
spectrum of the molecule in a linear-response measurement. In particular, we find that the phonon 
features are more pronounced for weak lead-molecule coupling. This way of measuring the excita- 
tion spectrum is less invasive than the more conventional current-voltage characteristic, which, by 
contrast, probes the system far from equilibrium. 
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PACS numbers: 73.23.Hk, 73.63.-b, 73.50.Lw, 85.65.+h 



I. INTRODUCTION 

The physics of electronic transport through single 
molecules has gained much interest in recent years, 
both experimentally and theoretically. Experiments with 
molecules ranging from H2 1 to DNA 2 have shown vari- 
ous interesting effects occurring in two-terminal 1-4 and 
three-terminal 5-9 molecular devices. The spectrum in- 
cludes Coulomb blockade, 7,9 negative differential resis- 
tance (NDR), 4 phonon influences, 6 ' 8 ' 9 and the Kondo 
effect. 7,9 The ultimate goal of these efforts is to realize the 
vision of molecular electronics. 10 While the possibilities 
of concrete applications remain to be seen, the question 
of transport in the molecular regime is of fundamental 
physical interest. 

Recently, there has been considerable theoretical 
effort to calculate current-voltage (IV) characteris- 
tics of single-molecule devices, and effects such as 
NDR, 11-13 influences of phonons and dissipation, 13-17 
and contact-geometry effects 18 ' 19 have been studied in- 
tensely. Roughly speaking, present approaches separate 
into two main directions: (1) Work based on a detailed 
modeling of the molecule and the contact region via den- 
sity functional theory (DFT). 18,19 For a recent discus- 
sion about the validity of equilibrium DFT for transport 
through single molecules cf. Ref. 20. (2) Approaches 
based on a parametric modeling of relevant molecular 
levels. 11-17 This type of approach has the advantage that 
it enables the investigation of additional degrees of free- 
dom including, e.g., mechanical, and magnetic degrees 
of freedom leading to phonon and spin dynamics. The 
present paper follows this second approach. 

In this paper, we investigate the thermopower of single- 
molecule devices. The thermopower is defined as the ra- 
tio of voltage V and an applied temperature difference 
AT under the condition that the current vanishes: 



S = — lhn -r-— 

AT^O AT 



(1) 
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For quantum dots the thermopower for pure sequen- 



tial tunneling has been investigated theoretically by 
Beenakker and Staring 21 and experimentally by Star- 
ing et al. 22 The cotunneling regime 23 and the crossover 
have been studied by Turek and Matveev. 24 In the case 
of a quantum dot strongly coupled to one lead, the 
thermopower of quantum dots has been investigated by 
Matveev and Andreev. 25 

Here, we extend these considerations to transport 
through single molecules, where experimental work 6,8,9 
indicates that phonons may play an important role. In 
our model, we consider transport through relevant elec- 
tronic orbitals and incorporate Coulomb interaction by a 
Hubbard-like term. These electronic orbitals are coupled 
to both oscillations of the molecule relative to the leads 
and internal vibrations. The coupling of the molecule 
to the leads is represented by a tunneling Hamiltonian. 
Employing the rate-equation approach valid for weak 
molecule-lead coupling, we compute the thermopower as 
a function of gate voltage, temperature, and electron- 
phonon coupling. 

We find that the thermopower contains information on 
the electronic and phononic excitations of the molecule. 
This way of measuring the molecular excitations in lin- 
ear response (cf. the AT — > limit in Eq. (1)) may have 
advantages over the more conventional TV characteris- 
tic. The latter necessarily involves nonequilibrium ef- 
fects, which are difficult to interpret. Moreover, a large 
applied voltage may affect symmetry and structure of the 
molecule itself. 

In leading order perturbation theory for the molecule- 
lead coupling, electrons tunnel from a lead onto the 
molecule or vice versa (sequential-tunneling contribu- 
tions). For pure sequential tunneling we find that the 
thermopower as a function of gate voltage develops a saw- 
tooth behavior in the low temperature limit with steps 
due to electronic and phononic excitations. Step sizes 
and their dependence on the electron-phonon coupling 
strength are analyzed. 

Moreover, we find that in a wide range of parameters 
so-called cotunneling contributions from next-leading or- 
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der perturbation theory for the molecule-lead coupling 
are important. 26 In this case, the electron only virtu- 
ally occupies molecular levels. We investigate the cotun- 
ncling contributions and the full crossover between the 
sequential-tunneling and the cotunneling regimes. We 
find that elastic cotunneling does not show any signifi- 
cant phonon structure, and discuss under which condi- 
tions sequential-tunneling phonon features are retained 
in the total thermopower. 

The outline of this paper is as follows: Section II in- 
troduces our model for a single-molecule device with me- 
chanical degrees of freedom. In section III we review the 
rate equations approach and address the issue of regu- 
larization of the cotunneling contributions. Our calcu- 
lations for the thermopower are described in section IV 
and the results arc presented in section V. We summa- 
rize our findings in section VI. Some calculational details 
including the results of the cotunneling regularization are 
relegated to appendices. 

II. MODEL 

We consider a three-terminal single-molecule device, 
consisting of a molecule weakly coupled to two metallic 
leads serving as source and drain electrode, respectively. 
The third electrode only influences the molecule by elec- 
trostatic interaction and serves as a gate electrode. In 
order to measure the thermopower of the device, it is 
necessary to control the electrostatic potentials and tem- 
peratures of the source and drain electrodes individually. 

The model we apply in order to investigate the ther- 
mopower has previously been used in analyses of IV 
characteristics, see, e.g., Refs. 14,15. The Hamiltonian 
can be divided into a part describing the electronic and 
phononic features of the molecule, a part modeling the 
leads, and a tunneling term that couples molecule and 
leads, H = H mo \ + H lcll ds + #mix, where 

H mo i =(e - eV g )n d + ^n d (n d - 1) 

+ \tku vih (tf + b)n d + hu vib (tfb + 1/2) (2) 

-Hleads = J]] J]] e P C apcr C a.pcr , (3) 
a=L,R p,er 

ff mi x= J2 E M z ) c U d - + h - c • ( 4 ) 

a=L,R p, a 

In the following, the common Fermi energy of the leads 
at vanishing bias voltage is chosen as the zero point of 
energy. For simplicity, we assume that only a single spin- 
degenerate orbital of the molecule with one-particle en- 
ergy £ contributes to the current. (A generalization to- 
wards more orbitals is not difficult.) For double occu- 
pancy of the molecule, Coulomb blockade is taken into 
account via the charging energy U . The operator d a 



(d^) annihilates (creates) an electron with spin projec- 
tion (7 on the molecule, n d — ^ CT dJ.rf CT denotes the 
molecule occupation-number operator. The whole sys- 
tem of molecular levels can be shifted by means of ap- 
plying a gate voltage V g . 

The leads are described as a non-interacting Fermi gas 
of electrons with a constant density of states. Here, 
c apa ( c apa) annihilates (creates) an electron in lead a 
(a = L, R) with momentum p and spin projection a. The 
potential V a and temperature T a of the left and right lead 
are taken into account through the probability distribu- 
tions for state occupation in the leads. It is assumed that 
relaxation in the leads is sufficiently fast so that at any 
time these distributions have the form of Fermi functions: 

f a (E) = (exp[(E + eV a )/k B T a ] + l)" 1 . (5) 

We distinguish two types of phonons, which we term vi- 
brations and oscillations: Vibrations are internal phonon 
modes of the molecule, for which the center of mass 
(CM) of the molecule is at rest, while oscillations in- 
volve movement of the molecule as a whole. Vibrational 
phonons are annihilated (created) by b (tf). For oscil- 
lations we use the momentum and position operators p z 
and z of the CM displacement. In the case of physisorbed 
molecules 31 the coupling to the leads is weak, so that 
these two phonon types typically involve different en- 
ergy scales: Vibrations, which are associated with strong 
intra-molccular bonds, will have considerably higher en- 
ergies than oscillations. The two phonon types also differ 
in the nature of coupling: Vibrations directly couple to 
the electric charge on the molecule, described by the term 
<~ n d {b^ + b), whereas the coupling for oscillations oc- 
curs through displacement-dependent tunneling matrix 
elements t^^z). 

(a) Oscillations. — Since t arises due to tunneling pro- 
cesses between the leads and the molecule, we assume an 
exponential fall-off of t with increasing distance between 
lead and molecule. For a symmetric molecule of length 
21 between two leads with a separation distance 2d, this 
yields 

th,R{z) = t exp[-(d - I ± z)/z ]. (6) 

The parameter zo fixes the length scale of the exponential 
fall-off of the electronic wave functions outside the leads 
and the molecule. 

(b) Vibrations. — For the vibrational electron-phonon 
coupling there exists a procedure which eliminates the 
coupling term by a canonical transformation of the 
Hamiltonian. 14,27 This yields a renormalization of the 
parameters e and [/, and of the lead-molecule coupling 
t(z) — > t(z) exp[— A(frt — b)]. Henceforth, for the sake of 
notational simplicity, we will refer to the renormalized 
parameters as e and U . (Alternatively, one can diagonal- 
ize the Hamiltonian for each occupation number n d and 
calculate Franck-Condon matrix elements.) 

In the following, we restrict ourselves to considering 
one phonon type at a time. Whenever the specific phonon 
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type is irrelevant we will skip the subscripts "vib" and 



"osc" 



III. 



RATE EQUATIONS AND TRANSITION 
RATES 



We consider the weak-coupling regime for the mole- 
cule-lead coupling. In this case the energy broadening 
7 of molecular levels due to i? m ; x is the smallest energy 
in the problem. In particular, we assume 7 <C fc^T, fkv, 
which allows for a perturbative treatment for H m { x . 

In the absence of coupling to the leads, the eigenstates 
of the molecule can be written as | na, q ) , where n de- 
notes the number of additional electrons on the molecule 
and q gives the number of excited phonons. The spin ori- 
entation a =|, I is only relevant for the singly-occupied 
molecule. Since we consider a spin-degenerate orbital 
and spin-independent tunneling matrix elements tL,n(z), 
cf. Eq. (2) and (4) , there exists a symmetry between the 
two states | 1 f,?} and 1 1 i,q). This allows for a nota- 
tionally more simple treatment without reference to spe- 
cific spin states by introducing appropriate spin factors 
into the transition rates, which account for the multiplic- 
ity of the n = 1 level, cf. Appendix A. 

The operator 7? m ; x introduces transitions between the 
eigenstates \n,q), for which the rates are calculated via 
Fermi's golden rule in the next subsection. Subsequently, 
these are used to formulate the rate equations and the 
expression for the steady-state current. 



A. Transition rates 

We abbreviate the total rate for a transition | n, q } — > 
\n',q') by W'^^f . These total transition rates can be 
written as a product of a Fermi factor f a or (1 — f a ), 
which gives the probability for the availability of elec- 
trons or holes at the appropriate energy in lead a, and a 
"bare transition rate" factor T^Z^qi calculated by Fermi's 
golden rule. We denote the energy of the molecule in the 
state I n, q } by 



E r ' 



n(s - eV g ) + Un(n - l)/2 + M<? + V 2 )- ( 7 ) 



Leading order perturbation theory yields sequential- 



n + !,</) 





FIG. 1: Sequential-tunneling processes, here schematically ex- 
emplified for tunneling between the molecule and the left lead. 



tunneling processes | n, q ) — > \n±l,q'), cf. Fig. 1, with 
the total rates 



W, 



w, 



n — >(n— 1) 



£ /. - K) r ^ 1} . (8) 

=L,R 

J2 [l-f a (E--E^ 



-iti — >(n— 1) 
q^q';a • 



a=L,R 



(9) 



The bare transition rates T are obtained by using Fermi's 
golden rule, 



■ q^q';a 



x |(n-l,g',e a |i7 mix |n,<7,0; 



M. 



n — >(n— 1) 



(10) 

2 



and analogously 



F 



n— >(n+l) 
q^q';a 



(n+1) 



K) 



M 



m(nti) 



q^q';a 



(11) 

Here, p a denotes the density of states in lead a, and 
s n^ m denies the spin factor, cf. Appendix A. In our 
calculations we assume Pl = PR = const. We note 
that due to our choice of symmetric tunneling matrix 

elements t^iz) = tii(—z), cf. Eq. (6), all rates r™^, n±1 ^ 
are de facto independent of the lead index o. The ma- 
trix elements M™^™^ 1 ^ can be expressed in terms of 
Laguerre polynomials and their detailed form depends 
on the phonon type, see Appendix B. 



\n,q) -» \n,q) 



\n,q) -» Kg') 





FIG. 2: Cotunneling processes, a) Elastic cotunneling, b) 
inelastic cotunneling 

Next-leading order perturbation theory generates co- 
tunneling processes with a virtual intermediate molecule 
state. 26 Cotunneling leaves the electron number on the 
molecule unchanged. By contrast, the phonon state can 
be changed, and in this case the process is called inelas- 
tic cotunneling. If initial and final state of the molecule 
are identical, one speaks of elastic cotunneling, cf. Fig. 2. 
In addition to the examples depicted in Fig. 2, where 
the electron is transferred from the left to the right lead, 
there are also contributions for the reverse process and 
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for processes which involve tunneling back and forth be- 
tween only one lead and the molecule. Thus, we abbre- 



viate the total cotunncling rates by W" 



q — >q'; a — >b 



where a 



and b denote cither L or R for the left and right lead. 

In order to obtain the total rates, one needs to sum 
over all lead energies e a weighted by Fermi functions for 
occupancy probabilities in the appropriate leads: 



(12) 



= J de a Pa (e a )f a (e a ) [l - f b (e b )} T" q ^, :a ^ b (e a ) , 



£ b = e a + - £™, 
vation. The bare cotunneling rates are given by 



where e = e a + E™ — E™, as required by energy conser- 



-pn— m 

a- 



>(0 



h 



Pb (e a -E" + E") 



j,n->(n+l) M (n+l)^n 
1V1 q^q":a 1V1 q"^q';b ^ 



M 



(13) 

2 



q^q":b 



M q»->q>;a 



e a + E% 



Ep 1 



E n q , 



e;,t 



Here, the first term in the sum represents cotunncling 
processes with an electron virtually tunneling from a lead 
into the molecule and out again. Conversely, the second 
term contains those cotunneling processes in which an 
electron virtually tunnels out of the molecule and another 
electron tunnels back in subsequently Again, the sym- 
metry of the tunneling matrix elements, th(z) = tji(-z), 



pn — m 

1 q^q';R- 



and FT*",. 



leads to r^.^ 

1 q^q': R^R- 

It has been pointed out previously in the literature that 
due to the singularities of the bare transition rates (13) 
for cotunneling, the integrals (12) for W diverge at finite 



temperatures. 
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Physically, this problem arises since 



we have assumed a well-defined energy for the interme- 
diate virtual state. In reality, due to higher order tun- 
neling effects, this intermediate state gains a finite width 
Ti, which leads to a regularization of the integrals. In 
the limit I\ — > one can derive a regularization scheme 
analogous to that presented by Turek and Matveev in 
Rcf. 24. This scheme consists of two steps: Firstly, intro- 
duce finite widths Ti of the intermediate states by adding 
an imaginary part in the energy denominators. Secondly, 
in order to avoid double counting, subtract terms scaling 
as 1/Tj, which correspond to sequential-tunneling con- 
tributions. In the r, — > limit, the remaining integrals 
can be written as Cauchy principal- value integrals. For- 
mally, this scheme can be extended to our case including 
phonons. However, this turns out to be a rather awkward 
procedure concerning the numerical treatment of the in- 
tegrals. Instead, we find it useful to transform the inte- 
grals so that a power series expansion of the expressions 
with respect to is feasible. After the subtraction of all 
1/Tj terms, the 1^ — > limit can be evaluated explicitly 
and we can express our results in terms of polygamma 
functions, cf. Appendix C. 



B. Rate equations 

In the weak-coupling regime and for nondegenerate 
eigenstates of the molecule, it can be shown that the 
off-diagonal elements of the reduced density matrix of 
the molecule p mo i are negligibly small. 14 ' 29 Therefore, the 
generalized master equations obtained via a density ma- 
trix approach reduce to simple rate equations. Writing 
Pq(t) = (n,q\ p m oi(i) \ n,q) for the probability that the 
molecule is in the state \ n,q) at time t, one obtains: 



dPl l /dt 



EE 

n'jtn q'^q 



Y q' q'^q 



pn\jjn^n 



(14) 



We consider the leading and next-leading order con- 
tributions to the transition rates derived in the previ- 
ous subsection. Generally speaking, sequential tunneling 
is the dominant process close to the Coulomb peaks in 
dl/dV, i.e., whenever alignment of a molecular level with 
the Fermi energies of the leads permits tunneling: 



k B T > mm \eV g 

m=0,l 



mU\ 



(15) 



Cotunncling plays the dominant role in the Coulomb val- 
leys and for sufficiently low temperatures, i.e., when 



ksT <C min \eV a — e — mU\ . 

m=0,l 



(16) 



Using the results for the transition rates, one obtains for 
the stationary case 

n = ^_pn _ [p(n-l) Tir(n-l)->n , p(n+l) w (n+l)^n 

Qt i Zj [ ?' V->9 q' q'^q 



q' 



— P n W™^(™ +1 ) _ p" 

i q^q' q q^q' 



+ \p n , w n ,^ n - p n w n_> ?l 

' / / L q' "q'^q 1 q ,r q^q'i 

q'^q 



[pn pcq pn 
r q r q Z^iqi q' 



(17) 



where we have included an additional term, which takes 
into account relaxation of the phonons in the relax- 
ation time approximation with relaxation time r. P| q = 
g-gfiw/fesT^ _ e -hw/k B Tj ig the equ iiib r hmi probability 
distribution for the phonons alone. The various cotunnel- 
ing contributions have been abbreviated by introducing 



W n 

q^q 



?= E E 

i=L,Rb=L,R 



q~*q ; a ~ 



(18) 



The stationary-state rate equations in conjunction with 
the normalization condition ^ n P™ = 1 form an in- 
homogeneous system of linear equations whose solution 
gives the stationary probability distribution P™ for a 
given voltage bias or temperature difference. 
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The stationary current is given by 



w 



ra— >(n+l) 
q^q';R 



-n — >(n— 1) 
q^q';R 



1 L yv q^q'-R^L ty q^q'; 



(19) 



L^R\ 



The first sum comprises all sequential-tunneling contri- 
butions, the second sum the cotunneling contributions. 

Calculations of IV characteristics based on this ap- 
proach show phonon steps in the IV curve for the sequen- 
tial tunneling regime and phonon steps in dl/dV for the 
cotunneling regime at low temperatures. 14,15 The relax- 
ation time approximation provides a means of analyzing 
situations between the two extremes of equilibrated and 
non-equilibrated phonons studied by Mitra et al. 14 



IV. THERMOPOWER 

The thermopower, Eq. (1), is calculated by considering 
the current through the molecule in the linear response 
regime, which is 

J(V, AT) = GV + G T AT + 0(V 2 , AT 2 , VAT), (20) 

where G denotes the conductance and Gt the thermal 
coefficient. Hence, the thermopower can be written as 



c _ = ZJL 

G ~ G S( i + G co 



Geo 
T 



(21) 



where sequential tunneling and cotunneling contributions 
have been separated. We investigate sequential tunneling 
and cotunneling contributions to the thermopower and 
obtain expressions valid in the full crossover regime by 
means of the regularization scheme (Appendix C). 

In order to obtain G and Gt, we expand the current 
(19) in the bias voltage V = Vl — Vr and the temperature 
difference AT = Tj, — Tr between the source and drain 
electrodes. Since V and AT are in principle infinitesi- 
mal, we can conveniently choose the right electrode to 
have zero potential and temperature T. Accordingly, the 
left electrode has potential V and temperature T + AT. 
When expanding the current, one has to expand both the 
probabilities P™ and the transition rates W™^} . 

We write the expansion for the transition rates and 
probabilities as 



rr q^q' 



w q^q> 



ATt^ 



q^q' 



and 



Pq = Pq + ©^AT + $™F + 



(22) 



(23) 



Here, P™ = 2* 1 '" exp(-P"//c B T)/Z denotes the 
grandcanonical probability distribution, and Z = 
J2 n q % Sl ' n exp(— P"//csT) is the corresponding partition 



function. (The additional factor of 2* 5 ™ 1 takes the spin- 
degeneracy of the level n = 1 into account.) The nor- 
malization condition for the deviations of P™ from its 
equilibrium value is = J2 n , q « = Sn, g ®q ■ 

By substituting the expansions (22), (23) into the rate 
equations (17) and retaining only terms linear in V and 
AT, one obtains a new set of rate equations for the devi- 
ations 0™ and <J>™, see Appendix D, Eq. (Dl). These, in 
conjunction with the normalization conditions for and 
<f>, represent an inhomogeneous system of linear equa- 
tions, whose solution yields the deviations 0™ and $™. 

Finally, the current (19) can be expanded in terms of 
AT and V: 



(AT0™ + V*J) 



w 



n— >(ra+l) 
^q';R 



— W 



i->(n-l) 
->g'; R 



n,q,q' 



J pn \ \T(+ n ~ >n — j.n— >n 

\ r q l^ 1 \ L q->q'; R^L L q^q': 



L->R) 



(24) 



+V(v. 



n — >n 

q^q'; R— *L 



v q^q'; L->r)] 



terms proportional , . . . . 



q 



to 0™ and 
of the cotunneling contributions 
they come with the coefficient 



Note that the 
in the expansion 
are absent since 

[wqZtq'-R-^L ~ w q^q'- l->r\i which vanishes due to the 
symmetry r^. a ^ 6 (e) = r^. 6 ^ a (e). 

(a) Sequential-tunneling contributions. — In Eq. (24) it 
was chosen to expand Ir in V and AT. Due to the 
steady-state property J = Ir = II, an expansion in II 
gives the same result and it turns out to be convenient 
to expand the expression I = (II + Ir)/2, which gives 
for the sequential-tunneling contributions 



n.q.q' 



-(n+l) 
*q';R 



-(n-1) 

*q'\ R 



+w 



♦ (n-1) 
q^q';L 



— W„ 



-(n+l) 
*q';L 



+ - 



+ 



n,q,q' 
n,q,q' 



pny 

q 



P"AT 



♦(n-1) 



— V„ 



>q>;L 



-(n+l) 
*q';L 



-(n+l) 
*q';L 



(25) 



Here, the first term remarkably vanishes due to the sym- 

„ ,. ; ~q^q™^L^- Therefore, one obtains the 
following sequential-tunneling contributions to the ther- 
mal coefficient Gt and the conductance G, 



metry w a ^, R 



r< s i 

\Jrp 



= \Y. 



pn 

q 

n,q,q> 
v pn 

q 

n,q,q> 



,n— >(n— 1) 
J q^q';L 



-t: 



-(n-1) 



-(n+l) 
*q';L 

-(n+l) 



(26) 



(27) 



We point out that the so-obtained conductance and ther- 
mal coefficient do not depend on the probability devia- 
tions 0^ and $™ any more. This is an important re- 
sult since it allows for an analytic expression of the ther- 
mopower not involving an explicit solution of the rate 
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equations, cf. Appendix D. (Expansions of Jj, and Ir 
alone lead to expressions for G sq and G^?, which do in- 
volve G™ and $™. We have also carried out calculations 
based on this approach by solving the rate equations for 
the probability deviations and find agreement with the 
results from Eqs. (26) and (27).) 

(b) Cotunneling contributions — The cotunneling con- 
tributions to thermal coefficient and conductance are 

Gt = E - t£?;^*), (28) 

n,q,q' 

G co = £ Tf(v£%. v n^n ^ (2g) 

n,q,q' 

In principle, any set of rate equations involving phononic 
excitations yields an infinite system of linear equations. 
In numerical calculations one makes use of the fact 



that transition matrix elements involving highly excited 
phonon states are typically very small. This allows for 
the introduction of a cutoff phonon number. 

We find that the linear response quantities G and Gt 
do not depend on the relaxation time r. Mathematically, 
this corresponds to the result that the conductance and 
thermal coefficient do not involve the probability devi- 
ations and 6™, cmp. Eqs. (26)-(29). The physical 
reason for this is the following: In the J — > limit the 
average time needed for one electron tunneling through 
the molecule becomes large compared to the relaxation 
time. Consequently, the initial state for any tunneling 
process corresponds to an equilibrium phonon state. 

By substituting back Eqs. (26)- (29) into Eq. (21), we 
arrive at the following analytical expression for the ther- 
mopower: 



T , P n 

£—in,q,q' q 



,n — y(n— 1) 
t q-,q';L 



,71— >(n+l) 



* L q^q';R^L ZjL q^q'-L^R 



L~in,q,q ! q 



V 



n — >(n— 1) 



v 



*(n+l) 

q^q'-L 



2v r : 



; R^L 



q—*q ; L — >H 



(30) 



This equation is our central result. It shows that even 
in the presence of phonons the thermopower can be ex- 
pressed analytically through the equilibrium probabil- 
ity distribution P™ and the expansion coefficients of the 
transition rates evalutated at vanishing source drain volt- 
age and temperature difference. In the following section 
the implications of Eq. (30) will be discussed. 

V. RESULTS 
A. Sequential tunneling 

We first consider the results for pure sequential tun- 
neling, postponing the discussion of the full thermopower 
due to both sequential and cotunneling to Sec. VB. We 
give numerical results for the thermopower and present 
analytic expressions for the limiting case U — ► oo and 
T — > below. Representative numerical results are 
shown in Fig. 3. 

The dominant feature of the sequential-tunneling ther- 
mopower S sq is a large step at the gate voltage V* = 
(e + U /2)/e for which the Fermi energy of the leads lies 
halfway in between the 1 1, ) and the | 2, } state. This 
situation is depicted in Fig. 4. At V g = V* the ther- 
mocurrent GtAT vanishes due to electron-hole symme- 
try. An increased (decreased) gate voltage lowers (raises) 
the molecular levels with respect to the Fermi level, 
and therefore current is dominated by electrons (holes) 
flowing from the left to the right lead. Consequently, 
Gt changes sign at V g = V*. Moreover, away from 



the Coulomb peaks, sequential tunneling can only occur 
through the tails of the lead Fermi distributions due to 
energy conservation. Therefore, the sequential-tunneling 
conductance and thermal coefficient fall off exponentially 
away from the Coulomb peak. At gate voltages close to 
V* , their behavior can be estimated by 

G sq ~ cxp[-e(y; - e)/k B T] cosh[e(V; - V g )/k B T], 

(31) 

Gy ~ exp[-e(y; - e)/k B T] sinh[e(y; - V g )/k B T). 

(32) 

Consequently, the sequential-tunneling thermopower 
S s<1 = Gy q /G sq behaves as - t&nh[e(V* - V g )/k B T] in 
the vicinity of V* , which develops a discontinuity in the 
limit T -> 0. 

In addition to the electronic step, the results show 
smaller phonon steps with a distance of fko between ad- 
jacent steps. To understand slope, temperature depen- 
dence, and phonon step sizes of the sequential-tunneling 
thermopower, we turn to the case U — > oo. 

In the U — > oo limit, electronic double occupation of 
the molecule is forbidden, and the sequential-tunneling 
thermopower S sq = G^?/G Bq reads 

ssq T, q , q > gTga - y "') _ Vg- e/e 

T.q,q'Pq\P°V^- Pl ^) ^ 

(33) 





FIG. 3: Thermopower S sq times temperature as a function 
of gate voltage. The parameters are chosen as e = 0, hu — 
0.11U. 



FIG. 5: Phonon step size Aq in units of htu/e for the step 
numbers Q — 1, 2, 3. 



V g and T S have units oi U/e. Large: Vibrations Qth gtep of T £ Sq [n thc T ^ q limit from Eq (34) j 



with A = 2. Inset: Oscillations with £0 = 5. (See text below 
Eq. (35) for a discussion of the choice of parameters.) The 
positions of the corresponding Coulomb peaks in 81/ dV are 
marked with arrows. In contrast, the main features of the 
thermopower occur between the Coulomb peaks. 



1 = 2 



n = 1 



T L = T + AT 



T„ = T 



FIG. 4: Electron-hole symmetry for the gate voltage Vg 
(e + U/2)/e. 



Thus, we find that the thermopower purely due to se- 
quential tunneling roughly scales like 1/T, which is in 
agreement with the quantum dot case. 21 In the low tem- 
perature limit, the thermopower as a function of gate 
voltage develops a characteristic sawtooth behavior. The 
slope of the linear pieces is found to be dS sq /dV g = 1/T. 
In the T — > limit, one obtains 



lim TS sq 



V g - e/e 



(34) 



- sgn(eVg - e) 



<\eVg-e\/tu*> 



e J2q<\eV g -e\/hu r 0^J 



where the last term generates the step features by adding 
up higher phonon contributions for increasing gate volt- 
ages. We can obtain the phononic step size Aq of the 



A Q =hoj/eTlZQ 
Q 



g=0 



|_<j=0 9=0 



(35) 



Here, Q counts the discontinuities of TS sq starting at 
e — eVc = with increasing gate voltage. Aq depends on 
the step number Q and on the coupling strength A or £ = 
Zo/Aosc for vibrations or oscillations, respectively. Here, 
Aosc = (fi/Mwosc) 1 / 2 is the harmonic-oscillator length for 
oscillations and M denotes the molecular mass. 

It is instructive to estimate typical values of the pa- 
rameters £ an d A for realistic systems. For £ we need 
to compare z to the oscillator length A osc . An order of 
magnitude estimate yields zq ~ h/ (nieW) 1 ^ 2 . Here, W is 
the work function of the metal leads and is of the order of 
several eV. On the other hand, for a typical experiment 6 
oscillations occur on an energy scale of 1-10 meV. This 
yields £0 = z /\ OS c 3> 1. In this case, displacements 
of the molecule's CM are small on the scale of zq, and 
therefore no significant shuttle effects can be expected. 

Next we consider the vibrational coupling parame- 
ter A. Let r be the normal coordinate deviation from 
the equilibrium value tq. An additional electron has 
to lowest order the effect of shifting the phonon po- 
tential curve by some distance Ar, so that the poten- 
tial energy is now w ^Mu>^ ih (r + nAr) 2 . Hence, the 
electron-phonon coupling term is of the order of magni- 
tude of Mu)l ih rAr = Ar/A v ibfiwvib(6 + b^) and therefore 
A ~ Ar/A v ib. Here, A v ;b = {h/ Mlj^) 1 / 2 is the harmonic 
oscillator length corresponding to vibrations. There is 
not a general rule for how A v ib and Ar compare so that 
A can in principle assume values both smaller and larger 
than 1. 

Due to the different behavior of the matrix elements 
for vibrations and oscillations, the phonon step size 
turns out to differ between those two cases as shown 
in Fig. 5. For vibrational phonons the electron-phonon 
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coupling becomes stronger for increasing A. In the case 
of electron-phonon coupling for oscillations, the coupling 
gets stronger for decreasing £ = Zo/\ osc (decreasing Zq 
at fixed A osc increases the position dependence of the 
hopping matrix elements t(z j). Thus, the plausible find- 
ing is that in both cases phonon step size increases with 
electron-phonon coupling strength. For oscillations the 
steps are rather small in the relevant regime of £o 3> 1- 
For vibrations they may be more pronounced. 



B. Results for the total thermopower 

The results for the thermopower discussed above arise 
from considering sequential-tunneling contributions only. 
However, if the Fermi levels are not aligned with a molec- 
ular level, sequential tunneling only occurs via electrons 
(or holes) in the tails of the Fermi distributions in the 
leads. In this case, the sequential-tunneling conductance 
G sq and thermal coefficient G^S are exponentially sup- 
pressed, and higher-order processes such as cotunncling 
may yield important contributions. Accordingly, sequen- 
tial tunneling dominates in proximity to the aligncd- 
levels configuration, and is suppressed most at the gate 
voltage V* at which the large electronic step in the 
sequential-tunneling thermopower occurs. In the latter 
range of gate voltages, cotunneling may give the domi- 
nant contributions to the thermopower. For this reason 
we have included the effect of cotunneling processes in 
the rate-equations approach, cf. Sec. III. At tempera- 
tures ksT < fad inelastic cotunneling can be neglected, 
and we find that the elastic cotunneling does not exhibit 
significant phonon structure. 

Figure 6 exemplifies the behavior of the thermopower 
including both sequential and cotunncling as a function of 
gate voltage for several temperatures. Whether the total 
thermopower S shows the sequential-tunneling phonon 
structure or whether it is mainly dominated by cotun- 
neling contributions without significant phonon features, 
strongly depends on the choice of parameters. Firstly, 
step-like features can only be expected if the sequential- 
tunneling part develops pronounced steps. As discussed 
above, this depends on phonon type, phonon-coupling 
strength, and temperature. Only for temperatures well 
below Hio/kB one can expect any features, as can be seen 
from the smoothening of the phonon steps for increasing 
temperature in Fig. 6. 

Secondly, temperature and the dimensionless coupling 
parameter, 

a = pt 2 Q /U, (36) 

which arises in the rate equations and roughly de- 
scribes the relative strength of cotunneling to sequen- 
tial tunneling, 32 determine where the crossover between 
the sequential-tunneling and cotunneling regimes occurs. 
For illustration, Fig. 7 shows the thermopower S as a 
function of gate voltage at fixed temperature for two dif- 




V 

g 



FIG. 6: Thermopower times temperature as a function of gate 
voltage for several temperatures. The parameter choices are: 
s = 0, hw — 0.11L/, A = 2 (vibrations), a — 1CP 9 (dimension- 
less coupling parameter defined in Eq. (36)). V g and TS are 
plotted in units of U/e. 





























total thermopower, a- 10 5 

total thermopower, a-10 9 

sequential-tunneling thermopower 





0.5 1 1.5 

V 

g 

FIG. 7: Thermopower as a function of gate voltage at fixed 
temperature. The parameter choices are: e — 0, fko = 0.11(7, 
A = 2 (vibrations), T = Q.QlU/ks- V g is plotted in units of 
U/e, S in units of fcs/e. 



ferent coupling parameters a as well as the corresponding 
sequential-tunneling result for comparison. 

The crossover between the sequential-tunneling and 
cotunneling regimes occurs in a rather small gate volt- 
age range, cf. Fig. 6 and 7, which allows one to identify 
crossover gate voltages V^ . 33 Our results show that the 
crossover points roughly scale as V*° ~ Tina -1 , which 
is in agreement with corresponding results for quantum 
dots, cf., e.g., Ref. 24. 

This can be understood based on the following es- 
timate of the crossover points V™. We assume that 
only a small gate-voltage region is dominated by se- 
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quential tunneling, V*° < V*. With eAV g — 
min TO= o > i \eV g — e — mU\ being the dominant activation 
energy for either electrons or holes, one can roughly esti- 
mate the sequential-tunneling conductance and thermal 
coefficient by an activated behavior dependence, 

G sq , G s T q ~ exp(-eAV g /k B T). (37) 

While sequential-tunneling contributions therefore fall off 
exponentially with AV g , cotunncling contributions only 
show a weak power-law dependence on the activation en- 
ergy eAVg, and temperature T. To lowest order they 
may be approximated by a constant, 

G co ,G c T °~ab. (38) 

Comparison of equations (37) and (38) yields as an esti- 
mate for the crossover gate voltage 

V™ » k B T(ln a- 1 - In b) /e. (39) 

In this approximation, V™ increases linearly with tem- 
perature and decreases logarithmically with lno; -1 . The 
number of phonon steps (if present in the sequential- 
tunneling contribution) is given by eV g ° /ftu). 

For the parameter choices of Fig. 7, we find that b 
assumes values so that the In b term in Eq. (39) can be 
neglected. This leads to the following estimates of the 
crossover gate voltage: V™ ~ O.lU/e for a = 1CP 5 , and 
Vg° w 0.2U/e for a = 10~ 9 , which is in good agreement 
with the crossovers observed in Fig. 7. 

We note that the thermopower attains rather small 
values in the cotunncling regime. This is plausible when 
reconsidering the effect of the electron-hole symmetry at 
eV g * = (e+U /2). Due to this symmetry, the thermopower 
must vanish at V* . In the case of sequential tunneling, 
the exponentially suppressed current G S ^AT shows the 
breaking of this symmetry for small gate voltage devia- 
tions from V* rather abruptly (leading to the large steps 
in the sequential-tunneling thermopower S sq as a func- 
tion of gate voltage). For cotunncling on the other hand, 
the thermal current is not exponentially suppressed, but 
roughly follows a power-like decrease with AV g . There- 
fore, breaking of the electron-hole symmetry is not as 
pronounced so that S remains small in the cotunncling 
region centered around V* . 

VI. SUMMARY 

Using a model for electronic transport through a sin- 
gle spin-degenerate molecular orbital, and taking into ac- 
count oscillational and vibrational phonons, we have cal- 
culated the thermopower of a single-molecule device in 
the regime of weak molecule-lead coupling. In contrast to 
IV measurements, the thermopower provides a means of 
extracting information about electronic and phononic ex- 
citations and the nature of the electron-phonon coupling 
in a linear response measurement. Therefore, it may have 



advantages over the more conventional IV characteris- 
tic, which necessarily involves noncquilibrium effects, and 
which, at large voltages, may even affect symmetry and 
structure of the molecule itself. 

We have found that sequential-tunneling contributions 
yield a characteristic sawtooth behavior of the ther- 
mopower as a function of gate voltage for low tem- 
peratures, which shows structure due to electronic and 
phononic excitations. It has been shown that due to 
the different nature of electron-phonon coupling for os- 
cillations and vibrations, characteristic differences in the 
phonon step size arise in the sequential-tunneling ther- 
mopower. Analytical expressions for the phonon step size 
have been derived in the limit of strong Coulomb block- 
ade, which show that, for realistic parameters, phonon 
steps can be expected to be more pronounced for vibra- 
tions than for oscillations. 

Away from the Coulomb peaks, cotunneling domi- 
nates over sequential tunneling. 26 By considering cotun- 
ncling contributions obtained by means of a regulariza- 
tion scheme, we have found that the (elastic) cotunncl- 
ing regime does not show significant structure due to 
phononic excitations. We have investigated the crossover 
regime and have given an estimate for the gate voltage 
at which the crossover takes place. It has been shown 
that the phononic structure exhibited by the sequential- 
tunneling contributions is retained in the thermopower if 
k B T <C hjo and ksThia^ 1 > huj, i.e., (1) the tempera- 
ture is low enough so that the phononic structure is not 
blurred out, and (2) the dimensionless coupling parame- 
ter a is so small that the crossover gate voltage is high 
enough to allow for at least one phonon feature in the 
thermopower. 
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APPENDIX A: SPIN FACTORS 

In our model we consider transport through one spin- 
degenerate molecular orbital. Including the phononic ex- 
citations, a basis of the corresponding Hilbert space is 
given by the states 

|0,g), |1T,«>, |U,«>, |2,g>. (Al) 

Due to the spin-degeneracy and the spin-independent 
tunneling matrix elements, one can make use of the re- 
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suiting symmetry and obtain a description without ex- 
plicit reference to spin indices. To see this, we write 
down the rate equations including spin indices, 

d 

pna 

dt q 



q>,*> 



(ri-l)cr' yy-(n-l)a'- 

g'^g 



>na j~,{n+l)(j' Ti/C n +l) (j/ ~^na 
+ q' W q'^q 



■Tna—>(n+l)a f 
q^q' 



pna yh-na^(n-l)a' 
g q^q' 



+ 



+ 



g'^g 



na TJina—^na pna TJina—*na 

q' —>q q q-^>q' 



(A2) 



q' .a'^a 



pna' 



jxrna — ma 



pna Tirna — mer 

g g^q' 



where it is implied that a = for n = 0, 2 and a =|, J. for 
n = 1. Here, PF denotes total transition rates between 
specific spin states. 

r(l±l)-l = 

*g' 



By defining P^ 1 = + P^ and 



wf^p^ 1 '* + the sequential-tunneling con- 

tributions to the rate equations can be cast into the form 
of Eq. (17). Accordingly, the sequential-tunneling spin 
factors are given by 



,(i±i)- 



1, 



2. 



(A3) 



For cotunneling transitions W g ™^ 00 and W 2 ^, 20 , one 
has to take into account that the virtual intermediate 
state is now spin-degencrate, which leads to a factor of 
4 due to the coherent sum, cf. Eq. (13). For cotunnel- 
ing transitions with n = 1, the intermediate state has no 
spin-degeneracy, but the rate equations contain an addi- 
tional spin- flip channel, cf. Eq. (A2), leading to a factor 
of 2. Therefore, by defining the cotunneling spin factors 

by 



s CW0 = s 2^2 



(A4) 



APPENDIX B: MATRIX ELEMENTS 

(a) Matrix elements for oscillations. — The matrix ele- 
ments M£^*7? Q for oscillations are given by: 



M 



M' 



92! 



t (q'\ cxp[-z/z ] \q) 

1/2 / \ 92-91 



(Bl) 



1 



e 4 «o Lg- 



f 1 ) 



-(nil) 
*q';R 



t (q'\ cxp[+z/z ] I q) 
1/2 



(B2) 



e*S L %-« (-^) 



where qi = min{q, q'}, q2 = max{g, q'}, and £0 = 
•ZoAosc- Finally, L"„(i) denotes the generalized Laguerre 
polynomial. Note that there is no dependence on whether 
the final electronic state is n + 1 or n — 1. 



(b) Matrix elements for vibrations- 
obtains, cf. Rcf. 14, 



For vibrations one 



„^(„-l) 
9^9 ; a 



to(q'\e- X ^-^\q)=to( qi \/q2\) 1/2 
x X^-^e-^l 2 L®-' 1 (A 2 ) 




for q' > q 
for </ < g, 
(B3) 



where again (71 = min{g, g'}, q 2 = max{g, g'}. Note that 
there is no dependence on the lead index a. The corre- 
sponding matrix clement for a transition n — > n + 1 can 
be obtained by using 



,,n->(n-l) _ „ n ->(n+l) 
q—>q';a q'^q:a 



(B4) 



the rate equations involving spin, Eq. (A2), can be trans- 
formed into the set of rate equations in Eq. (17). 



See Sees. II and V A for definitions and a discussion of 
the parameters A OS c, zq and A. 



J 



APPENDIX C: REGULARIZATION SCHEME 



We sketch the regularization scheme that we apply in order to extract the cotunneling contributions from the 
diverging next-leading order perturbation theory. The integrals w, t, and v for cotunneling, cf. Eq. (22), contain 
factors of Fermi functions and their derivatives, energy factors and the bare transition rates. Our starting point is to 
introduce a finite width for the intermediate state. The bare transition rates then assume the following form: 



pro— »n 

1 9— >q' ; a^b 



fc=r 

£ 



e-E k +ir k ' e-E 



+ 



B k \ 

E' k +ir k J 







2' 


A k | 2 1 


B k 




\€-E k +iT k \ 


e-E' k +ir k 





+ 2 Re 



9 fc<9 



£-Bfc+ir fc e-E Q -ir q e-E' k +ir k e-E' g -iT q 
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Accordingly, we need to consider several types of integrals. In the following, we always evaluate the integrals in the 
r — > limit. The notation "— C(l/r)" indicates that terms proportional to 1/r have been subtracted before carrying 
out the limit. These terms correspond to sequential-tunneling contributions, cf. Ref. 15. 

1 = J de.f(e-E 1 )f(e-E 2 ) 7 ^ rrT7=7 ^ 

= (i7r + n B (-Bi -^ 2 )[-V(l/2 + i/3[^i -ei]/27r)+V(l/2-i/3[^i -e 2 ]/27r)] (CI) 

ei — £2 L 

+ n B (E 2 - E,) [-0(1/2 + i(3[E 2 - ei]/27r) + V(l/2 - \f3[E 2 - e 2 ]/27r)] | 
J = J de f(e - E{)f{e - E 2 ) _ * + p2 + " - 0(l/r)" 



^ {n B (^i - £2) Im (1/2 + i/3[£ - E 1 ]/2ir) + n B (E 2 - Pa) Im V> (1) (1/2 + i/3[E- P 2 ]/2tt)} 



(C2) 



Here, ip^(x) denotes the polygamma function of order n, cf. Ref. 30, and Ub(x) — (exp(x/fcsT) — 1) _1 is the Bose 
function. By evaluating I (J) in the limit E 2 — > 00, one obtains an expression for the integral I' (J') with only one 
Fermi factor. 

K = J def'(e E,)f{e i? 2 )_L__L_ = - dI/dEl (C3) 

L = J def'(e - Bx)/(c - ^ 2+r 2 = - 9J / 9 ^ ( C4 ) 

M = I de{e-E 1 )f(e-E 1 )f(e-E 2 ) e _^_ iTie _ e2 \ iT2 =G+{e 2 -E 1 )K (C5) 

M' = Jde(e- E 1 )f{e E,) [1 - /(e - E 2 )} £ _ £i ^ iri 6 _ 6a 1 +ira - /? 8/801' M (C6) 

iV = I de (e -£i)/'(e -£i)/(e - E 2 ) - _ ^ - p2 + « - P(l/r)" = Re G + {E - E X )L (C7) 

' = 1 & (e - ^/'(e - E,) [1 - /(e - £ 2 )] _ ^ + p2 + " - 0(l/r)» = 8/8(3 J' - N (C8) 

G = [ def(e - E x )f{e - E 2 ) l —- (C9) 

J e - ei - ll 

= - E 2 )^ l \l/2 + if3(E 1 - ei)/27r) + n'^ - P 2 ) [^(1/2 + i/3(Pi - ei)/27r) - ^(1/2 + i(3(E 2 - e 1 )/2 7 r)} 

Ztt 



N 



APPENDIX D: RATE EQUATIONS FOR 6™, EXPANSION COEFFICIENTS 

The rate equations for the deviations 0™ and <I>™ from the equilibrium probability have the following form: 







E 



+ ATt 



i->(n-l) 



Vv. 



'q'^q ) ' " <J' V— - "(J'^iJ 1 ' "<?'^<J ^ " <? V— - "g->g' 1 ' "q-^q' 

(n-l)^ _l MT«»-1 _I_ T/^-l^.t™- 1 )^™ _L ^AT«™+1 _I_ T/rf>"+l'U„( Il + 1 )- > ™ 



-(Are; i + ^)( % "4" +1 ) + 



w 



n-»(n-l) 
tj->q>' 



E 

<?'#<? 

i _ 



P^(ATt^rr q + Vv^q) - P?(*TtZ:$ + Vv^>) + (AT6™, + 



(AT0™ + V^)w 



n — > 7 
q) w q^qi 



Are™ + - p° q ^(ato;;, + v*j) 



(Dl) 
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By expanding (8) and (9) one obtains the following expressions for the expansion coefficients of the sequential-tunneling 
rates = £„ W^ 1 ': 



rw(n+l) _ n->(n+l) _ , / ™+i ™\ r 
q—>q'; L - W g^g';R -JRK^q' ~ h q ) 1 

1 - - 



? n\ r n^(n+l) 
9^9' 



n— *(n— 1) ro— >(n— 1) 



,n— >(n+l) 
l g^q';L 



g-^g'\ R 



n->(ra-l) 
L q^q':L 

ra— >(n+l) 
n— >(n— 1) 



tn-l\ pi- >(«-!) 
9' /J 9^9' ' 

- E%)/Tf' R (E% +1 E%) r^ n+1) , 

n-l\ r n— >(n-l) 



= +(££ - Ey)/Tf R {E n q - is*- 1 ) r; 1 

= +e/^ +1 - J B ? ")r^" +1) , 



Similarly, the expansion of the cotunneling rates (12) yields 

<^;«^6 = P / [l-/fl(e + ^-^)]r^ io ^ 6 (e) 



4.71— >n 



de- 



e + E? - E", 



q -fR{z)f R {t + E n q - E^)T n q ^ ;R ^ L (e) 



v 



,n— >(n+l) _ 
l q^q':R ~ 
,n — ► (n— 1) 
t q^q';R 
ji— >(n+l) 

9^9'; -R 

n— >(n— 1) 
' J q^q';R 



o. 

o, 
o, 

0. 



q—^q'i R^R 



u q->q'; L—>R 



u q->q'; R-^L 



v g-^g'\ L^L 



V 



T 

n—>n 



q^q'lR^R 







-ep | def' R (e) [l - / fl ( e + - E n q ,)] T^, ;L ^ R (e) 
ep J de f R (e)f R (e + E n q E n ql ) T£%. R ^ L (e) 

ep J de (f R (e)f' R (e + E n q - ££) f' R (e) [l - f R (e + E n q - E%)] ) T n q ^, L ^ L (e) 



(D2) 
(D3) 
(D4) 
(D5) 
(D6) 
(D7) 

(D8) 
(D9) 
(D10) 



de ( * +E \ E *' f R (e)f R (e + E n q - E$) - j,f' R (e) [l - f R (e + E n q E n q ,)] ) r££. L ^(e) (Dll) 
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